library(knitr)
library(dplyr)
library(sf)
library(mapsf)
library(RColorBrewer)
library(leaflet)
library(htmlwidgets)
library(htmltools)
library(ggplot2)
library(plotly)
library(potential)
library(car)Régression spatiale
INTRODUCTION
Objectifs
Cette séance va permettre :
de revoir plus en détail le concept de potentiel en utilisant le package R potential qui en facilite le calcul et la cartographie.
de revoir les modèles de régression multiple et de régression Poisson en les appliquant à la prédiction du nombre de clubs de sport présents ou absent d’une commune.
de mélanger des variables explicatives de type endogène (caractéristiques internes de la commune) et exogènes (situation de la commune par rapport aux espaces envrionnants)
Chargement des packages
L’installation est a priori la même que dans les sessions précédentes. On reprend juste pour mémoire la liste, à l’intention de ceux qui n’auraient pas suivi la séance précédente. On y ajoute un package pour le calcul de potentiel (potential) et un autre pour l’analyse des résultats de régression (car)
DONNEES
On reprend les trois jeux de données précédents en y ajoutant les deux bases de données sur le nombre de licences et le nombre de clubs présents dans chaque commune.
## Contour des communes
com<-readRDS("ParisPC/mapcom_parisPC.RDS")
## Population sur grille de 200 m
gri<-readRDS("ParisPC/gridpop_parisPC.RDS")
## Equipements localisés en latitude longitude
equ<-readRDS("ParisPC/equip_ParisPC.RDS")
## Nombre de licences
lic <- readRDS("ParisPC/lic2018_ParisPC.RDS")
## Nombre de clubs
clu <-readRDS("ParisPC/clu2018_ParisPC.RDS")Tableau de variables endogènes
On extrait de chaque fichier les informations relatives au sport retenu pour l’analyse et on centralise les résultats dans le fichier communal. On choisit de travailler sur l’exemple du football en prenant comme code de fédération 111 et comme code d’équipement 2802.
# Zone d'étude
tabcom <- st_drop_geometry(com)
# Extraction des clubs
myfede<-c("111")
myclu<- clu %>% filter(code_federation %in% myfede) %>%
select(insee_com, nbclu =clubs_sportifs_2018)
# Extraction des licences du sport choisi
myfede<-c("111")
mylic<- lic %>% filter(code_fed %in% myfede) %>%
group_by(insee_com) %>%
summarise(nblic=sum(nblic))
# Extraction du nombre total de licences
lictot<- lic %>% group_by(insee_com) %>%
summarise(lictot=sum(nblic))
# Extraction des équipements
typequ <- c("2802")
myequ <- equ %>% st_drop_geometry() %>%
filter(typ_code %in% typequ) %>%
group_by(insee_com) %>%
summarise(nbequ=n())
# Variables socio-eco
soceco<- gri %>% st_drop_geometry() %>%
group_by(insee_com) %>%
summarise(Ind = sum(Ind),
Men = sum(Men),
Men_pauv = sum(Men_pauv),
Ind_snv = sum(Ind_snv)) %>%
mutate(pop = Ind,
men = Men,
txpauv = 100*Men_pauv/Men,
revhab = Ind_snv/Ind) %>%
select(insee_com, pop,men, txpauv, revhab)
# Assemblage des tableaux
tabcom<- tabcom %>% left_join(soceco) %>%
mutate(denpop = pop/sup) %>%
left_join(lictot) %>%
left_join(mylic) %>%
mutate(pctlic = 100*nblic/lictot) %>%
left_join(myclu) %>%
left_join(myequ) %>%
arrange(insee_com)
# Remplacement des NA par des 0
tabcom$nbclu[is.na(tabcom$nbclu)]<-0
tabcom$nblic[is.na(tabcom$nblic)]<-0
tabcom$nbequ[is.na(tabcom$nbequ)]<-0
# données + fonds de carte
mapcom <-left_join(com, tabcom)CARTOGRAPHIE
Prenons l’exemple d’une cartographie du nombre de licenciés (stock) et de la part du total des licences (taux).
Avec mapsf
Normalement vous connaissez toutes les fonctions utilisez. Bien noter le changement de projection de la carte en crs = 2154.
mapcom<-mapcom %>% st_transform(2154)
mapdep<-mapcom %>% group_by(dept) %>%
summarise()
# Choix des classes
mybreaks<-c(0, 1,2,4,8,16,32,100)
# Choix de la palette
mypal<-brewer.pal(7, "RdYlBu")
# Carte de taux (choroplethe)
mf_map(mapcom, type="choro",
var="pctlic",
breaks= mybreaks,
pal=mypal,
border="white",
lwd=0.5,
leg_title = "% total licence",
leg_val_rnd = 0,
leg_pos = "topright")
# Ajout des départements
mf_map(mapdep, type="base",
col=NA,
border="black",
lwd=1,
add=T )
# Carte de stock (proportionelle)
mf_map(mapcom, typ = "prop",
var = "nblic",
inches = 0.05,
col="gray50",
leg_title = "nombre de licences",
leg_pos = "topleft")
# Cadre, titre, ..
mf_layout(title = "Distribution des licences de golf dans Paris PC",
credits = "Source : INSEE et Min. des Sports",
frame = T,
arrow=F,
scale=T)
Avec leaflet
Programme plus compliqué … mais qui permet d’aboutir au même résultat avec une carte interactive. Noter que la projection n’est pas modifiée et demeure crs = 4326.
map<-mapcom %>% select(dept,insee_com, nom_com, nblic, pctlic ) %>% st_transform(4326)
map$lng<-st_coordinates(st_centroid(map))[,1]
map$lat<-st_coordinates(st_centroid(map))[,2]
mapdep<-mapcom %>% group_by(dept) %>%
summarise() %>%
st_transform(4326)
# Choix de la variable
myvar <-map$pctlic
# Choix des classes
mybreaks<-c(0, 1,2,4,8,16,32,100)
# Choix de la palette (c'est une fonction !)
mypal <- colorBin('RdYlBu',
myvar,
bins=mybreaks)
# Calcul du diamètre des cercles
myradius <-8*sqrt(map$nblic/max(map$nblic, na.rm=T))
# Préparation des popups
mypopups <- lapply(seq(nrow(map)), function(i) {
paste0( paste("Commune : ",map$nom_com[i]), '<br>',
paste("Code INSEE : " ,map$insee_com[i]), '<br>',
paste("Nb. de licences : " ,map$nblic[i]), '<br>',
paste("% total licence :", round(map$pctlic[i],1))
)
})
mypopups<-lapply(mypopups, htmltools::HTML)
map <- leaflet() %>%
addProviderTiles('Esri.WorldTopoMap') %>%
# Réalisation de la carte choroplèthe
addPolygons(data = map,
fillColor = ~mypal(pctlic),
fillOpacity = 0.5,
color = "white",
popup = mypopups,
weight = 1,
highlightOptions = highlightOptions(weight = 3, color = 'green')) %>%
# Ajout de la carte des départements
addPolygons(data = mapdep,
fill = FALSE,
color = "black",
weight = 2) %>%
# Ajout de la carte de stocks
addCircleMarkers(data=map,
lat = ~lat,
lng = ~lng,
radius = myradius,
stroke = FALSE,
label = ~nblic,
fillColor = "gray50",
fillOpacity = 0.5)%>%
# Ajout de la légende
addLegend(data = map,
pal = mypal,
title = "% licences",
values = ~pctlic,
position = 'topright')
mapREGRESSION LINEAIRE
On va ici un modèle de régression n’utilisant que des variables endogènes (internes aux communes). On essaye de modéliser Y = % de licences pour le sport considéré
Choix des variables
Au vu de la carte précédente, on peut penser à plusieurs variables explicatives et faire des hypothèses sur le résultat attendu
- X1 : densité de population (relation positive car il faut de la place pour installer un terrain de football)
- X2 : revenu moyen par habitant (relation négative car les riches pratiquent plutôt des sports d’élites ce qui n’est pas le cas du football)
- X3 : % de ménages pauvres (relation positive car le football est considéré à tort ou à raison comme un outil de promotion sociale)
On crée un tableau avec Y, X1, X2, X3
don<-tabcom %>% select(insee_com, nom_com, Y=pctlic, X1=denpop, X2=revhab,X3=txpauv)
head(don) insee_com nom_com Y X1 X2 X3
1 75101 PARIS-1ER-ARRONDISSEMENT 1.703085 8391.209 38433.52 13.23255
2 75102 PARIS-2E-ARRONDISSEMENT 5.564024 21350.505 34985.14 15.00129
3 75103 PARIS-3E-ARRONDISSEMENT 3.400690 30778.632 36195.03 13.99706
4 75104 PARIS-4E-ARRONDISSEMENT 4.714971 15908.750 36304.64 13.55867
5 75105 PARIS-5E-ARRONDISSEMENT 3.075031 19597.619 37720.86 13.18044
6 75106 PARIS-6E-ARRONDISSEMENT 4.521818 16633.721 43607.03 12.49411Analyse des corrélations
On analyse la matrice de corrélation :
cor(don[,3:6]) Y X1 X2 X3
Y 1.0000000 -0.35258137 -0.8156021 0.72652758
X1 -0.3525814 1.00000000 0.2526739 0.04524385
X2 -0.8156021 0.25267387 1.0000000 -0.72002742
X3 0.7265276 0.04524385 -0.7200274 1.00000000On vérifie la forme des relations et on teste leur significativité
- Densité de population
scatterplot(don$X1,don$Y)
cor.test(don$X1,don$Y)
Pearson's product-moment correlation
data: don$X1 and don$Y
t = -4.474, df = 141, p-value = 0.00001569
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.4884600 -0.2000084
sample estimates:
cor
-0.3525814 - Richesse par habitant
scatterplot(don$X2,don$Y)
cor.test(don$X2,don$Y)
Pearson's product-moment correlation
data: don$X2 and don$Y
t = -16.738, df = 141, p-value < 0.00000000000000022
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.8640692 -0.7521516
sample estimates:
cor
-0.8156021 - % ménages pauvres
scatterplot(don$X3,don$Y)
cor.test(don$X3,don$Y)
Pearson's product-moment correlation
data: don$X3 and don$Y
t = 12.555, df = 141, p-value < 0.00000000000000022
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.6385290 0.7957734
sample estimates:
cor
0.7265276 Toutes les variables affichent des corrélations linéaires significatives ! Mais on note que la forme des relations n’est pas forcément linéaire. On y reviendra …
Modélé additif
On teste le modèle additif suivant
\(Y = a_0+a_1.X_1+a_2.X_2+a_3.X_3+ \epsilon\)
modlin<-lm(Y ~ X1 + X2 + X3, data=don)
summary(modlin)
Call:
lm(formula = Y ~ X1 + X2 + X3, data = don)
Residuals:
Min 1Q Median 3Q Max
-7.7420 -1.9405 -0.2735 1.5522 12.0714
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 19.59836228 2.06344477 9.498 < 0.0000000000000002 ***
X1 -0.00020419 0.00003642 -5.607 0.000000107254 ***
X2 -0.00039364 0.00005677 -6.934 0.000000000141 ***
X3 0.35444421 0.05503644 6.440 0.000000001810 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.117 on 139 degrees of freedom
Multiple R-squared: 0.7598, Adjusted R-squared: 0.7546
F-statistic: 146.6 on 3 and 139 DF, p-value: < 0.00000000000000022anova(modlin)Analysis of Variance Table
Response: Y
Df Sum Sq Mean Sq F value Pr(>F)
X1 1 699.1 699.1 71.940 0.00000000000002938 ***
X2 1 3170.6 3170.6 326.279 < 0.00000000000000022 ***
X3 1 403.0 403.0 41.476 0.00000000180962499 ***
Residuals 139 1350.7 9.7
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- Conclusion : toutes les variables sont significatives et affichent les signes attendus. L’analyse de variance montre cependant que c’est la variable X2 (revenu par habitant) qui est la plus déterminante. La qualité de l’ajustement est élevée (r2 = 76%)
Modèle multiplicatif
Essayons maintenant un modèle multiplicatif de la forme
\(log(Y) = a_0+a_1.X_1+a_2.X_2+a_3.X_3+ \epsilon\)
qui correspond à une forme multiplicative de l’effet des variables puisque l’on a :
\(Y = exp(a_0+a_1.X_1+a_2.X_2+a_3.X_3+ \epsilon)\)
et donc :
\(Y = exp(a_0).exp(a_1.X_1).exp(a_2.X_2).exp(a_3.X_3)\)
modexp<-lm(log(Y) ~ X1 + X2 + X3, data=don)
summary(modexp)
Call:
lm(formula = log(Y) ~ X1 + X2 + X3, data = don)
Residuals:
Min 1Q Median 3Q Max
-1.2206 -0.1266 0.0045 0.1464 0.5606
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.817787144 0.157533985 24.235 < 0.0000000000000002 ***
X1 -0.000012995 0.000002781 -4.673 0.00000692 ***
X2 -0.000053438 0.000004334 -12.330 < 0.0000000000000002 ***
X3 0.007414592 0.004201765 1.765 0.0798 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.238 on 139 degrees of freedom
Multiple R-squared: 0.7949, Adjusted R-squared: 0.7905
F-statistic: 179.6 on 3 and 139 DF, p-value: < 0.00000000000000022- Conclusion : Pas de différences dans le signe des coefficients et leur significativité. Mais la qualité de l’ajustement est plus élevée (r2 = 79%) ce qui signifie que le modèle multiplicatif décrit plus fidèlement les effets constatés.
POTENTIELS
On va maintenant utiliser le package potential pour créer des variables de type exogène mesurant la distribution des indicateurs non pas dans la commune mais dans le voisinage de celle-ci.
Potentiel d’équipement
# Extraction des équipements et projection 2154
typequ <- c("2802")
myequ <- equ %>%filter(typ_code %in% typequ) %>%
mutate(nb = 1) %>%
st_transform(2154)
# projection des communes en 2154
mapcom <- mapcom %>% st_transform(2154)
# Distance communes-équipement
dist<-create_matrix(myequ,mapcom)
# calcul du potentiel
mapcom$pot_equ_2000<-potential(x=myequ, # Ressources
y=mapcom, # Population
d=dist, # Distance population x Ressources
var = "nb", # Quantité de ressource
fun="e", # famille exponentielle
span = 2000, # distance ou f(Dij) = 0.5
beta=2) # Exposant de la distance
# Cartographie du résultat
mf_map(mapcom, type="choro",var="pot_equ_2000")
mf_map(myequ, type="base",add=T,col="red")
- Commentaire : La carte montre que les arrondissements du centre de Paris qui ne disposent pas de terrain de football ont cependant un potentiel d’acès à ceux-ci dans un voisinage gaussien de 2 km.
Potentiel de licences
# Extraction des licences projection 2154
maplic<-mapcom%>% select(nblic) %>%st_transform(2154)
# projection des communes en 2154
mapcom <- mapcom %>% st_transform(2154)
# Distance communes-équipement
dist<-create_matrix(maplic,mapcom)
# calcul du potentiel
mapcom$pot_lic_2000<-potential(x=maplic,y=mapcom,d=dist,var = "nblic", fun="e",span = 2000,beta=2)
# Cartographie du résultat
mf_map(mapcom, type="choro",var="pot_lic_2000")
mf_map(maplic, type="prop",var="nblic",col="red", inches=0.05)
- Commentaire : La carte montre que le potentiel de sportifs ayant une licence de football est maximale dans les arrondissements périphériques de Paris et dans la banlieue Nord.
REGRESSION DE POISSON
On va construire un modèle de régression de Poisson pour prévoir le potentiel de clubs du potentiel d’équipement et de son potentiel de licence.
Choix des variables
Le tableau comporte uniquement des stocks et des potentiels
Y : nombre de clubs X1a : nombre d’équipements dans la commune X1b : potentiel d’équipement dans un voisinage gaussien de 2.5km X2a : nombre de joueurs ayant une licence dans la commune X2b : potentiel de joueurs ayant une licence dans un voisinage gaussien de 2.5km
On crée un tableau avec Y, X1a, X1b X2a,X2b
don<-mapcom %>% select(insee_com, nom_com,
Y=nbclu,
X1a = nbequ,
X1b = pot_equ_2000,
X2a = nblic,
X2b = pot_lic_2000) %>%
st_drop_geometry()
head(don) insee_com nom_com Y X1a X1b X2a X2b
1 75119 PARIS-19E-ARRONDISSEMENT 17 3 37.262588 2162 22122.189
2 75106 PARIS-6E-ARRONDISSEMENT 5 0 6.116627 400 14617.246
3 93071 SEVRAN 3 13 47.018154 1190 10356.941
4 94056 PERIGNY 2 1 5.746092 84 814.633
5 94041 IVRY-SUR-SEINE 5 10 38.133234 1305 17336.098
6 75112 PARIS-12E-ARRONDISSEMENT 14 29 65.061731 1393 28111.675Analyse des corrélations
On analyse la matrice de corrélation :
cor(don[,3:7]) Y X1a X1b X2a X2b
Y 1.0000000 0.4315624 0.2416476 0.7516748 0.5794876
X1a 0.4315624 1.0000000 0.6279578 0.5791498 0.2387389
X1b 0.2416476 0.6279578 1.0000000 0.5980269 0.4025822
X2a 0.7516748 0.5791498 0.5980269 1.0000000 0.5992941
X2b 0.5794876 0.2387389 0.4025822 0.5992941 1.0000000On vérifie la forme des relations et on teste leur significativité
- Equipement de la commune
scatterplot(don$X1a,don$Y)
cor.test(don$X1a,don$Y)
Pearson's product-moment correlation
data: don$X1a and don$Y
t = 5.6808, df = 141, p-value = 0.00000007394
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.2878018 0.5563023
sample estimates:
cor
0.4315624 - Potentiel d’équipement dans le voisinage de la commune
scatterplot(don$X1b,don$Y)
cor.test(don$X1b,don$Y)
Pearson's product-moment correlation
data: don$X1b and don$Y
t = 2.957, df = 141, p-value = 0.003643
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.08070004 0.39031390
sample estimates:
cor
0.2416476 - Nombre de licences dans la commune
scatterplot(don$X2a,don$Y)
cor.test(don$X2a,don$Y)
Pearson's product-moment correlation
data: don$X2a and don$Y
t = 13.533, df = 141, p-value < 0.00000000000000022
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.6702226 0.8152346
sample estimates:
cor
0.7516748 - Potentiel de licences dans le voisinage de la commune
scatterplot(don$X2b,don$Y)
cor.test(don$X2b,don$Y)
Pearson's product-moment correlation
data: don$X2b and don$Y
t = 8.4432, df = 141, p-value = 0.00000000000003385
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.4590001 0.6790441
sample estimates:
cor
0.5794876 Modèle
On teste le modèle suivant :
\(Y = exp(\alpha+\beta_{1a}.X_{1a}+\beta_{1b}.X_{1b}+\beta_{2a}.X_{2a}+\beta_{2b}.X_{2b} +\epsilon)\)
On emploie une régression de Poisson car le nombre de clubs est une variable quantitative discrète composée d’entier. On suppose qu’elle obéit à une loi de Poisson et on vérifie si c’est exact en comparant la moyenne et l’écart-type de Y.
mean(don$Y)[1] 3.79021sd(don$Y)[1] 3.711143Des tests plus précis pourraient être utilisés mais dans le cas présent on voit que les deux valeurs sont très proche et qu’il est acceptable de considérer que le nombre de club obéit bien à une loi de Poisson.
modpoi<-glm(Y ~ X1a + X1b + X2a + X2b, data=don, family="poisson")
summary(modpoi)
Call:
glm(formula = Y ~ X1a + X1b + X2a + X2b, family = "poisson",
data = don)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.25253750 0.16302149 1.549 0.121356
X1a 0.03330562 0.00952547 3.496 0.000471 ***
X1b -0.02185387 0.00406491 -5.376 0.0000000761 ***
X2a 0.00086511 0.00008965 9.650 < 0.0000000000000002 ***
X2b 0.00005579 0.00001095 5.096 0.0000003464 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 417.02 on 142 degrees of freedom
Residual deviance: 121.44 on 138 degrees of freedom
AIC: 534.46
Number of Fisher Scoring iterations: 5Anova(modpoi,type="III")Analysis of Deviance Table (Type III tests)
Response: Y
LR Chisq Df Pr(>Chisq)
X1a 11.676 1 0.0006331 ***
X1b 29.632 1 0.00000005223 ***
X2a 94.633 1 < 0.00000000000000022 ***
X2b 25.540 1 0.00000043329 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Discussion
Les quatre facteurs conditionnent bien l’apparition des clubs de football dans les communes et sont tous significatifs. Le nombre de terrain de football dans la commune a un effet positif sur l’apparition d’un ou plusieurs club dans la commune ce qui parait logique. Par contre la quantité de terrains de football dans les communes voisines a un effet négatif ce qui peut se comprendre comme un effet de concurrence. Quant au nombre de personnes détenant une licence, il a un effet positif à la fois dans la commune et dans les communes voisines.
La régression de Poisson n’offre pas directement de qualité d’ajustement comme dans le cas de la régression par la méthode des moindres carrés ordinaires. On peut néanmoins calculer un pseudo-R2 dit de Mc Fadden en calculant la part de déviance expliquée c’est-à-dire en effectuant le calcul suivant :
\(R^2_{McFadden} = \frac{NullDeviance-ResidualDeviance}{NullDeviance}\)
(417.02-121.44)/(417.02)[1] 0.7087909Dans notre exemple on trouve un \(R^2\) de Mc Fadden d’environ 71% ce qui est très satisfaisant et montre que la connaissance de la localisation des terrains et des licences permet dans une large mesure de prédire le nombre de clubs présents dans une commune.
PROLONGEMENTS
Pour en savoir plus sur les modèles de régression, le plus simple est de partir du petit billet introduction au GLM de Claire Della Vedova qui explique bien la différence entre le modèle linéaire classique et les modèles de régression logistique ou de régression de Poisson.
Ensuite … il faudra s’attaquer à un bon manuel de statistiques. Nous recommandons celui de Daniel J. Denis Univariate, Bivariate Multivariate Statistics using R qui est très pédagogique et fournit les programmes d’application en R ou en Python.